Setup R notebook Environment
knitr::opts_chunk$set(echo = TRUE)
Sys.setlocale("LC_CTYPE", "en_US.UTF-8")
[1] ""
packages = c("sp","sf","tidyverse","tmap","rjson", "rgdal", "geojsonio", "GISTools",'spatialEco', "raster","dplyr", "spatstat",
'maptools', 'leaflet')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}
Shapefile Object with Spatial data Import data for past 12 years
Requires translation of chinese character independetly (https://stat.ethz.ch/pipermail/r-sig-geo/2009-March/005323.html)
Taiwan follows EPSG 3826 (http://spatialreference.org/ref/epsg/twd97-tm2-zone-121/)
taiwan_ts_map_sp <- readOGR(dsn = "data/TAIWAN_TOWNSHIP", layer = "TOWN_MOI_1071226", stringsAsFactors=TRUE)
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\terence.tan.2015\Desktop\dangy\data\TAIWAN_TOWNSHIP", layer: "TOWN_MOI_1071226"
with 368 features
It has 7 fields
#taiwan <- readShapePoly("data/taiwan_data/")
taiwan <- readOGR(dsn = "data/taiwan_data", layer = "COUNTY_MOI_1071226", stringsAsFactors=TRUE)
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\terence.tan.2015\Desktop\dangy\data\taiwan_data", layer: "COUNTY_MOI_1071226"
with 22 features
It has 4 fields
#taiwan@proj4string<- CRS( "+init=epsg:3826 +proj=longlat +ellps=WGS84 +no_defs")
#taiwan.union <- aggregate(taiwan)
df_dengue <- jsonlite::fromJSON("data/dengue_case.json") %>%
filter(as.Date(Onset_day) >= "1999-01-01" & as.Date(Onset_day)< "2000-01-01")
# Check for NA values for coordinates
sum(is.na(df_dengue$Minimum_statistical_area_center_point_X))
[1] 0
sum(is.na(df_dengue$Minimum_statistical_area_center_point_Y))
[1] 0
df_dengue<- df_dengue[!is.na(df_dengue$Minimum_statistical_area_center_point_X),]
df_dengue<- df_dengue[!(df_dengue$Minimum_statistical_area_center_point_X == 'None'),]
# Type conversion
df_dengue[, c(10,11,19,23,24)] <- sapply(df_dengue[, c(10,11,19,23,24)], as.numeric)
NAs introduced by coercionNAs introduced by coercion
# Transform into SF object
sf_dengue <- st_as_sf(df_dengue,
coords = c("Minimum_statistical_area_center_point_X",
"Minimum_statistical_area_center_point_Y"),
crs = "+init=epsg:3826 +proj=longlat +ellps=WGS84 +no_defs",na.fail=FALSE)
sf_dengue <- na.omit(sf_dengue)
sp_dengue <- as(sf_dengue, 'Spatial')
sp_dengue <- as(sp_dengue, 'SpatialPoints')
tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(taiwan)+
tm_borders(col = "grey40",alpha=0.5)+
# prevent zoom too much
tm_view(alpha = 1, set.zoom.limits = c(7,21))+
tm_layout(basemaps = c('OpenStreetMap'))+
tm_shape(sp_dengue)+
tm_dots(col ="red",size = 0.01)
As of version 2.0, basemaps and basemaps.alpha have to be called from tm_view
tmap_mode("plot")
tmap mode set to plotting
dengue_ppp <- as(sp_dengue, "ppp")
summary(dengue_ppp)
dengue_ppp <- as(sp_dengue, "ppp")
summary(dengue_ppp)
Planar point pattern: 699 points
Average intensity 148.3408 points per square unit
*Pattern contains duplicated points*
Coordinates are given to 6 decimal places
Window: rectangle = [120.15043, 121.79933] x [22.322371, 25.180104] units
Window area = 4.71212 square units
plot(dengue_ppp)
qt <- quadrat.test(dengueTW_ppp,
nx = 20, ny = 15)
plot(dengueTW_ppp)
plot(qt, add = TRUE, cex =.1)
tw_owin <- as(taiwan_ts_map_sp, "owin")
plot(tw_owin)
summary(tw_owin)
Window: polygonal boundary
1007 separate polygons (no holes)
enclosing rectangle: [118.13797, 124.56115] x [21.8956, 26.385278] units
Window area = 3.26423 square units
Fraction of frame area: 0.113
dengueTW_ppp = dengue_ppp[tw_owin]
dengueTW_ppp <- unique(dengueTW_ppp)
summary(dengueTW_ppp)
Planar point pattern: 646 points
Average intensity 197.9024 points per square unit
Coordinates are given to 6 decimal places
Window: polygonal boundary
1007 separate polygons (no holes)
enclosing rectangle: [118.13797, 124.56115] x [21.8956, 26.385278] units
Window area = 3.26423 square units
Fraction of frame area: 0.113
#plot(dengueTW_ppp)
#dengueTW_ppp.km <- rescale(dengueTW_ppp, 1000, "km")
kde_taiwan_bw <- density(dengueTW_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
#plot(kde_taiwan_bw)
plot(kde_taiwan_bw)